In [12]:
import numpy as np
from vivarium import Vivarium
from spatio_flux import PROCESS_DICT, TYPES_DICT, build_path, render_path

Make a Vivarium¶

view available types and processes

In [13]:
# TODO -- make a spatio-flux specific Vivarium with the the core preloaded
vi = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)
In [14]:
# view the available types
vi.get_types()
Out[14]:
['',
 'length/time^2',
 'mass/length*time',
 'mass^0_5/length^0_5*time',
 'length^2*mass/temperature*time^2',
 'length^3*mass/current^2*time^4',
 'length^2*mass/current*time^3',
 'length^3/mass*time^2',
 '/time',
 'tree',
 'number',
 'kinetics',
 'substance/length^3',
 'mass/current*time^2',
 'integer',
 'length^2*mass/current*time^2',
 'length^2*mass/time',
 'length^2*mass/current^2*time^3',
 '/substance',
 'current*time/mass',
 'length^3/time',
 'schema',
 'mass^0_5/length^1_5',
 'length^2*mass/time^2',
 'printing_unit',
 'current',
 'particle',
 'interval',
 '/length',
 'current*time',
 'positive_float',
 'length^1_5*mass^0_5/time',
 'length^0_5*mass^0_5/time',
 'length^4*mass/time^3',
 'length*mass/current^2*time^2',
 'string',
 'maybe',
 'enum',
 'current^2*time^3/length^2*mass',
 'printing_unit/length',
 'length',
 'current^2*time^4/length^2*mass',
 'length*temperature',
 'mass/length*time^2',
 'current*time/substance',
 'current*length^2',
 'length*mass/time^2',
 'edge',
 'current*length*time',
 'protocol',
 'path',
 'map',
 'luminosity',
 'boolean',
 'float',
 'length/time',
 'length*mass/current*time^3',
 'length^2/time',
 'substance/time',
 'any',
 'current*length^2*time',
 'length^2/time^2',
 'union',
 'length^2',
 'substrate_role',
 'mass/time^3',
 '/temperature*time',
 'length^3',
 'length^1_5*mass^0_5/time^2',
 'current^2*time^4/length^3*mass',
 'mass/temperature^4*time^3',
 'length^2*mass/current^2*time^2',
 'substance',
 'luminosity/length^2',
 'time/length',
 'mass/length',
 'current*time^2/length^2*mass',
 'time^2/length',
 'bounds',
 '/printing_unit',
 'wires',
 'array',
 'tuple',
 'boundary_side',
 'process',
 'mass',
 'length/mass',
 'quote',
 'composite',
 'temperature',
 'mass/time^2',
 'step',
 'reaction',
 'list',
 'length^0_5*mass^0_5',
 'length*time/mass',
 'length^2*mass/time^3',
 'time',
 'length^2*mass/substance*temperature*time^2',
 'emitter_mode',
 'mass/length^3',
 'length^4*mass/current*time^3']
In [15]:
# view the available processes
vi.get_processes()
Out[15]:
['json-emitter',
 'console-emitter',
 'DiffusionAdvection',
 'ram-emitter',
 'Particles',
 'MinimalParticle',
 'composite',
 'DynamicFBA']

dFBA¶

Dynamic Flux Balance Analysis (dFBA) extends traditional Flux Balance Analysis (FBA) to model the dynamic behavior of metabolic networks over time, allowing for the simulation of growth and substrate consumption in a changing environment.

In [16]:
# inspect the config schema for the 'increase' process
vi.process_config('DynamicFBA')
Out[16]:
{'model_file': 'string',
 'kinetic_params': 'map[tuple[float,float]]',
 'substrate_update_reactions': 'map[string]',
 'biomass_identifier': 'string',
 'bounds': 'map[bounds]'}
In [17]:
vi.core.default(vi.process_config('DynamicFBA')) # TODO use default to
Out[17]:
{'model_file': '',
 'kinetic_params': {},
 'substrate_update_reactions': {},
 'biomass_identifier': '',
 'bounds': {}}
In [18]:
dfba_config  = vi.process_config('DynamicFBA', dataclass=True)  # TODO -- make this work

Error finding process DynamicFBA: 'type_parameters'
In [19]:
# TODO - enable get_dataclass to work with the new process
# dfba_config  = v.get_dataclass('DynamicFBA')
dfba_config = {
    "model_file": "textbook",
    "kinetic_params": {
        "glucose": (0.5, 1),
        "acetate": (0.5, 2)},
    "substrate_update_reactions": {
        "glucose": "EX_glc__D_e",
        "acetate": "EX_ac_e"},
    "biomass_identifier": "biomass",
    "bounds": {
        "EX_o2_e": {"lower": -2, "upper": None},
        "ATPM": {"lower": 1, "upper": 1}
    }
}
In [20]:
# make the vivarium
v = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)

# add a dynamic FBA process called 'dFBA'
v.add_process(name="dFBA",
              process_id="DynamicFBA",
              config=dfba_config,
              )
v.diagram(dpi='70')
Out[20]:
No description has been provided for this image
In [21]:
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]

for mol_id in mol_ids:
    v.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand()
    )

v.connect_process(
    name="dFBA",
    inputs={
            "substrates": {
                mol_id: build_path(path, mol_id)
                for mol_id in mol_ids}
        },
    outputs={
            "substrates": {
                mol_id: build_path(path, mol_id)
                for mol_id in mol_ids}
        }
)

# add an emitter to save results
v.add_emitter()

# TODO -- show_values does not work
v.diagram(dpi='70', show_values=True)
Out[21]:
No description has been provided for this image
In [22]:
v.set_value(path = ['fields', 'glucose'], value=10)
v.set_value(path = ['fields', 'biomass'], value=0.1)
field = v.get_value(['fields'])
print(field)

{'glucose': 10, 'acetate': 0.012004913732089673, 'biomass': 0.1}
In [23]:
v.save(filename='dFBA_t0')

Saved file: out/dFBA_t0.json
In [24]:
# run the simulation for 10 time units
v.run(interval=60)
In [25]:
v.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=[
        [
            '/fields/glucose',
            '/fields/acetate',
            '/fields/biomass'
        ]
    ]
)
Out[25]:
No description has been provided for this image

Spatial dFBA¶

In [26]:
mol_ids = ["glucose", "acetate", "biomass"]
path=["fields"]
rows = 2
columns = 1

# make the vivarium
v2 = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)
for mol_id in mol_ids:
    v2.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand(rows, columns)
    )

# add a dynamic FBA process at every location
for i in range(rows):
    for j in range(columns):
        dfba_name = f"dFBA[{i},{j}]"

        # add a process for this location
        v2.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config,
        )
        v2.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}},
            outputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}}
        )

# add an emitter to save results
v2.add_emitter()

v2.diagram(dpi='70')
Out[26]:
No description has been provided for this image
In [27]:
# change some initial values
v2.set_value(path = ['fields', 'glucose', 0, 0], value=10)
v2.set_value(path = ['fields', 'biomass', 0, 0], value=0.1)
field = v2.get_value(['fields'])
print(field)

{'glucose': array([[10.        ],
       [ 0.49791182]]), 'acetate': array([[0.55230278],
       [0.88987915]]), 'biomass': array([[0.1      ],
       [0.9479142]])}
In [28]:
v2.run(60)
In [29]:
v2.get_timeseries(as_dataframe=True)
Out[29]:
/global_time /fields/glucose /fields/acetate /fields/biomass
0 0.0 [[10.0], [0.4979118174000552]] [[0.552302778281881], [0.8898791463284117]] [[0.1], [0.9479141961152532]]
1 1.0 [[9.904761904761905], [0.024946499325416316]] [[0.5644108304779702], [0.14649946583879703]] [[0.1079432568708625], [1.0034766298215592]]
2 2.0 [[9.80200585245463], [0.0]] [[0.5773867179541111], [0.0]] [[0.11651530697082768], [1.0128449439227383]]
3 3.0 [[9.691145527078628], [0.0]] [[0.5912815487035986], [0.0]] [[0.12576552148631978], [1.0128449439227383]]
4 4.0 [[9.571550338512791], [0.0]] [[0.6061469258667948], [0.0]] [[0.13574706721093277], [1.0128449439227383]]
... ... ... ... ...
56 56.0 [[0.0], [0.0]] [[0.0], [0.0]] [[0.97895496621264], [1.0128449439227383]]
57 57.0 [[0.0], [0.0]] [[0.0], [0.0]] [[0.97895496621264], [1.0128449439227383]]
58 58.0 [[0.0], [0.0]] [[0.0], [0.0]] [[0.97895496621264], [1.0128449439227383]]
59 59.0 [[0.0], [0.0]] [[0.0], [0.0]] [[0.97895496621264], [1.0128449439227383]]
60 60.0 [[0.0], [0.0]] [[0.0], [0.0]] [[0.97895496621264], [1.0128449439227383]]

61 rows × 4 columns

In [30]:
# get a list of all the paths so they can be plotted together in a single graph
all_paths = []
for i in range(rows):
    for j in range(columns):
        # get the paths for this location
        location_path = []
        for mol_id in mol_ids:
            this_path = build_path(path, mol_id, i, j)
            rendered_path = render_path(this_path)
            location_path.append(rendered_path)
        all_paths.append(location_path)
# print(all_paths)
In [31]:
v2.plot_timeseries(
    subplot_size=(8, 3),
    combined_vars=all_paths
)
Out[31]:
No description has been provided for this image
In [32]:
v2.show_video()

Diffusion/Advection¶

This approach models the physical processes of diffusion and advection in two dimensions, providing a way to simulate how substances spread and are transported across a spatial domain, essential for understanding patterns of concentration over time and space.

In [33]:
vi.core.default(vi.process_config('DiffusionAdvection'))
Out[33]:
{'n_bins': (0, 0),
 'bounds': (0.0, 0.0),
 'default_diffusion_rate': 0.1,
 'default_diffusion_dt': 0.1,
 'diffusion_coeffs': {},
 'advection_coeffs': {}}
In [34]:
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
    'glucose',
    'acetate',
    'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
    'biomass': (0, -0.1)
}
path = ['fields']

v3 = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)

for mol_id in mol_ids:
    v3.add_object(
        name=mol_id,
        path=path,
        value=np.random.rand(n_bins[0], n_bins[1])
    )

v3.add_process(name='diffusion_advection',
               process_id='DiffusionAdvection',
               config={
                   'n_bins': n_bins,
                   'bounds': bounds,
                   'default_diffusion_rate': diffusion_rate,
                   'default_diffusion_dt': diffusion_dt,
                   # 'diffusion_coeffs': diffusion_coeffs_all,
                   'advection_coeffs': advection_coeffs,
               },
               inputs={'fields': ['fields']},
               outputs={'fields': ['fields']}
               )

# add an emitter to save results
v3.add_emitter()

v3.diagram(dpi='70')
Out[34]:
No description has been provided for this image
In [35]:
v3.run(60)
In [36]:
v3.plot_snapshots(n_snapshots=5)
No description has been provided for this image
In [37]:
v3.show_video()

COMETS (Computation Of Microbial Ecosystems in Time and Space)¶

COMETS combines dynamic FBA with spatially resolved physical processes (like diffusion and advection) to simulate the growth, metabolism, and interaction of microbial communities within a structured two-dimensional environment, capturing both biological and physical complexities.

In [38]:
bounds = (10.0, 10.0)
n_bins = (10, 10)
mol_ids = [
    'glucose',
    'acetate',
    'biomass'
]
diffusion_rate = 1e-1
diffusion_dt = 1e-1
advection_coeffs = {
    'biomass': (0, -0.1)
}
path = ['fields']

v4 = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)

# create the molecular fields
# for mol_id in mol_ids:
v4.add_object(
    name='glucose',
    path=path,
    value=np.random.rand(n_bins[0], n_bins[1])
)
v4.add_object(
    name='biomass',
    path=path,
    value=np.ones((n_bins[0], n_bins[1])) * 0.01
)
v4.add_object(
    name='acetate',
    path=path,
    value=np.zeros((n_bins[0], n_bins[1]))
)

# add a diffusion/advection process
v4.add_process(name='diffusion_advection',
               process_id='DiffusionAdvection',
               config={
                   'n_bins': n_bins,
                   'bounds': bounds,
                   'default_diffusion_rate': diffusion_rate,
                   'default_diffusion_dt': diffusion_dt,
                   'advection_coeffs': advection_coeffs},
               inputs={'fields': path},
               outputs={'fields': path}
               )

# add a dynamic FBA process at every location
for i in range(n_bins[0]):
    for j in range(n_bins[1]):
        dfba_name = f"dFBA[{i},{j}]"

        # add a process for this location
        v4.add_process(
            name=dfba_name,
            process_id="DynamicFBA",
            config=dfba_config,
        )
        v4.connect_process(
            name=dfba_name,
            inputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}},
            outputs={"substrates": {
                        mol_id: build_path(path, mol_id, i, j)
                        for mol_id in mol_ids}}
        )

# add an emitter to save results
v4.add_emitter()

v4.diagram(dpi='70',
           remove_nodes=[f"/dFBA[{i},{j}]" for i in range(n_bins[0]-1) for j in range(n_bins[1])]
           )
Out[38]:
No description has been provided for this image
In [39]:
v4.run(60)
In [40]:
v4.plot_snapshots(n_snapshots=5)
No description has been provided for this image
In [41]:
v4.show_video()
In [42]:
v4.plot_timeseries(
    subplot_size=(8, 3),
    query=[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0',
    ],
    combined_vars=[[
        '/fields/glucose/0/0',
        '/fields/acetate/0/0',
        '/fields/biomass/0/0',
    ]]
)
Out[42]:
No description has been provided for this image

Particles¶

In [1]:
# TODO -- remove this from final
import numpy as np
from vivarium import Vivarium
from spatio_flux import PROCESS_DICT, TYPES_DICT, build_path, render_path
In [2]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (20, 40)  # Number of bins in the x and y directions

v5 = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)

v5.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.4,
        'boundary_to_add': ['top']
    },
)
v5.connect_process(
    name='particle_movement',
    inputs={
        'fields': ['fields'],
        'particles': ['particles']
    },
    outputs={
        'fields': ['fields'],
        'particles': ['particles']
    }
)

v5.initialize_process(
    path='particle_movement',
    config={'n_particles': 2}
)

v5.add_emitter()
v5.diagram(dpi='70')
Out[2]:
No description has been provided for this image
In [3]:
v5.run(100)
v5_results = v5.get_results()
In [4]:
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif

# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
    v5_results,
    skip_frames=3,
    bounds=bounds
)

Saving GIF to species_distribution_with_particles.gif
In [5]:
v5.diagram(dpi='70')
Out[5]:
No description has been provided for this image
In [6]:
v5.save(filename='v5_post_run', outdir='out')

Saved file: out/v5_post_run.json

Minimal particle process¶

In [7]:
# load the previous document
# TODO -- make this work
# v6 = Vivarium(document='out/v5_post_run.json',
#               processes=PROCESS_DICT,
#               types=TYPES_DICT)
# v6.set_value('global_time', 100.0)
#
# # get rid of all the particles
# v6.set_value('particles', {})
#
# # add a glucose field
# v6.add_object(
#     name='glucose',
#     path=['fields'],
#     value=np.random.rand(n_bins[0], n_bins[1])
# )
#
# # add a single particle
# v6.initialize_process(
#     path='particle_movement',
#     config={'n_particles': 1}
# )
# v6.diagram(dpi='70')
In [11]:
bounds = (10.0, 20.0)  # Bounds of the environment
n_bins = (20, 40)  # Number of bins in the x and y directions

v6 = Vivarium(processes=PROCESS_DICT, types=TYPES_DICT)

# add a glucose field
v6.add_object(
    name='glucose',
    path=['fields'],
    value=np.random.rand(n_bins[0], n_bins[1])
)

v6.add_process(
    name='particle_movement',
    process_id='Particles',
    config={
        'n_bins': n_bins,
        'bounds': bounds,
        'diffusion_rate': 0.1,
        'advection_rate': (0, -0.1),
        'add_probability': 0.4,
        'boundary_to_add': ['top']
    },
)
v6.connect_process(
    name='particle_movement',
    inputs={
        'fields': ['fields'],
        'particles': ['particles']
    },
    outputs={
        'fields': ['fields'],
        'particles': ['particles']
    }
)

v6.initialize_process(
    path='particle_movement',
    config={'n_particles': 1}
)
In [12]:
v6.run(10)
v6_results = v6.get_results()
In [13]:
from spatio_flux.viz.plot import plot_species_distributions_with_particles_to_gif

# TODO -- integrate this method with vivarium
plot_species_distributions_with_particles_to_gif(
    v6_results,
    skip_frames=2,
    bounds=bounds
)

Saving GIF to species_distribution_with_particles.gif